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Abstract 

We have carried out extensive computer simulations of one-dimensional mod- 
els related to the low noise (solid-on-solid) non-equilibrium interface of a two 
dimensional anchored Toom model with unbiased and biased noise. For the 
unbiased case the computed fluctuations of the interface in this limit provide 
new numerical evidence for the logarithmic correction to the subnormal La 
variance which was predicted by the dynamic renormalization group calcu- 
lations on the modified Edwards- Wilkinson equation. In the biased case the 
simulations are in close quantitative agreement with the predictions of the 

2 

Collective Variable Approximation (CVA), which gives the same Lz behavior 
of the variance as the KPZ equation. 
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I. INTRODUCTION 



The nature of the interface separating two equilibrium phases, or more generally any two 
distinct bulk states of matter, is a problem of continuing interest. While there is in most 
cases some fuzziness in the transition region, giving rise to an intrinsic structure, the width 
of this is usually much smaller than (and therefore separable from) the fluctuations in the 
location of the interface. These fluctuations are fairly well understood both microscopically 
and macroscopically in (simple) equilibrium systems PJ , but much less is known about them 
in the larger context of nonequilibrium situations M. 

An important advance in the latter case was the introduction of equations of the Edwards- 
Wilkinson H and Kardar-Parisi- Zhang Q type to describe fluctuations in dynamically mov- 
ing interfaces. Solutions of these equations appear to describe in a quantitative manner the 
macroscopic fluctuations of a large class of such interfaces @. It is probably fair to say, 
however, that there are few microscopic systems (even simple ones) for which one can argue 
a priori, with mathematical rigor, that their behavior will be described by these equations. 

Several years ago two of us (JLL and ES), together with B. Derrida and H. Spohn [|5]||, 
introduced and studied the behavior of the low-noise-limit sharp interface separating the + 
and — phases in a simple model nonequilibrium system — the 2D Toom cellular automaton. 
Our semi-infinite interface was situated in the third quadrant of the square lattice and was 
anchored at the origin. Its precise location could be specified by a spin configuration on the 
positive semi-infinite one-dimensional integer lattice Z+. Fluctuations of the location of the 
interface at a distance L from the origin are then directly transcribed into fluctuations of 
the magnetization Ml = Y^i=i a % m t ne stationary state of the one- dimensional model. 

The dynamics of the semi-infinite ID spin model are as follows: starting with some initial 
configuration (containing a mixture of + and — spins) we pick a site % with rate 1 if Oi = — 1, 
and with rate A < 1 if o% — 1, and exchange <Tj with the spin <jj, where j is the first site to 
the right of % such that cxj = — o\. Remarkably, the model inherits the unidirectional nature 
of information flow in the Toom model and the stationary state for the first L spins can be 
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obtained exactly (no finite size effects) by restricting the dynamics to a system with L spins, 
with the provision that if site i is in the last block of spins, i.e., if <7j = a i+ i = . . . = o~l, then 
the spin Oi is flipped, i.e., changed to —Oi (with rate 1 or A). The model can be classified 
into two cases, the unbiased case with A = 1, and the biased case with A 7^ 1. 

Quantities of interest in the model include the average of the magnetization (Ml) and its 
variation Vl = ((Me, — {Ml)) 2 ) as functions of the system size L. By symmetry, (Ml) = 
in the unbiased case, and a simple computation, using the well justified assumption that 
spins far from the origin are statistically independent, shows that in general (Ml) — fiL 
with fi = The fluctuations have not been calculated exactly. Simulations reported in 

@ gave Vl ~ L v with v = 0.53 in the unbiased case and v = 0.57 for A = 1/4, but in fact 
approximate treatments discussed there suggested that the true exponents are v = 1/2 in 
the unbiased case, possibly with logarithmic corrections, and v = 2/3 in the biased case. 

Two approximation methods were introduced in ||. The first, the collective variable 
approximation (CVA), gives Vl° v ~ ^3/2L^ 2 in the unbiased case and V^ v — CL 2 ! 3 
(with C a computable constant) in the biased case. The second approach was based on 
a description of the interface by a non-linear stochastic diffusion equation of the EW |3| 
and KPZ [|J type. Analysis of this equation predicts how fluctuations of a (doubly) infinite 
uniform interface will grow in time; in the semi-infinite problem, the fixed boundary at 
the left of the system and the positive velocity of excitations convert this growth to a 
corresponding Independence of Vl- This approach predicts Vl ~ L 2 ^ 3 in the biased case 
(here 2/3 is the usual KPZ exponent). In the unbiased case the prediction depends on the 
growth of excitations in a modified KPZ equation with a third order, but no second order, 
nonlinearity; a subsequent renormalization group analysis of this problem ( |7|], see also ||) 
leads to the prediction Vl ~ L 1 / 2 log^ 4 (L/L )- 

The simulation results of did not include sufficiently large systems to distinguish the 
logarithmic behavior in the unbiased case from behavior L u , with v a power slightly higher 
than 1/2, or to verify the power 2/3 in the biased case. Simulations of J7| and [§] support 
the renormalization group conclusions for a uniform system but, since the passage from the 
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uniform to the semi-infinite system is somewhat heuristic, provide only indirect guidance for 
the latter. In the present note we describe new simulations on our original model for L < 2 19 
both for A = 1 and A = 1/4. We also did simulations on modified models and models with 
different bias values. Our results on the unbiased case give support to the predictions of J7| . 
Our results for the biased case seem to suggest that the CVA is a very good approximation. 

For computational efficiency, we used a type of "multi-spin coding" , similar to the tech- 
nique applied for simulating the repton model in P,[10|. The computationally intensive part 
of the code is written in bit operations like AND (A), exclusive OR (ffi) and NOT (->). This 
allows one to run 32 or 64 independent simulations in parallel by applying bit operations to 
four-byte or eight-byte integers of which each bit corresponds to an independent simulation. 

We introduce a coding variable rjj = (1 + &j)/2. To implement the algorithm, we first 
choose (as described below) a marker m, whose nonzero bits select the simulations in which 
updating will take place, then pick a random site j and flip the selected spins at this site: 

r)'j — m © rjj (1) 

Next we walk along the chain, starting at k = j + 1 and incrementing k by one at each step. 
In each simulation we flip the spin a k if it differs from <jj and if so we set the corresponding 
bit in the marker to zero: 

d = Vk © Vj (2) 
rj' k = rj k © (m A d) (3) 
m = m A -<d (4) 

We continue walking until either our marker m equals or we have reached the end of 
the chain. To insure that simulations running in different bits do not become identical, we 
initialize m as m = r, where r is a random bit-sequence with each bit taking the value 1 
with probability 1/2, and then repeat with a different starting site j and with m = ->r. This 
generates a large degree of independence between simulations corresponding to different bits. 
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To obtain the mask m we use the random number generator marsag [11]], which has 
good randomness properties for all the bits; for the site selection we use the random number 



generator ranmar [11], which has a long sequence and good spectral properties. 



In the absence of a bias (A = 1), the thermalization time to required was determined to 
be t = L 2 /8 spin exchanges. The autocorrelation time of the magnetization in the steady 
state is much smaller, growing as L 3 / 2 , as found also in ||, although the thermalization 
time suggests that some correlations must grow as L 2 . The presence of a bias increases the 
required thermalization time; for A = 1/4, to ~ L 2 /4 and for A = 1/8, to ~ L 2 /2. For 
A = 1/2 we took again to ~ L 2 /A as thermalization, although we could have taken a bit 
less. One run typically starts with a random spin configuration, which is evolved over time 
(n + l)to- The magnetization in each simulation (corresponding to a bit) is measured every L 
moves; for each interval of length t a separate histogram of the magnetization M, including 
all bit simulations, is constructed. The first histogram is discarded (thermalization) and from 
the remaining n histograms we obtain moments of the magnetization. Reported statistical 
errors in these moments are one standard deviation errors based on the assumption that 
these n histograms are statistically independent. Note that even if simulations in different 
bits are correlated, we still obtain a reliable error estimation. 



II. THE UNBIASED CASE 

In the absence of bias (A = 1) we now have data for sizes L < 2 19 , presented in Table 
I. Two analytic descriptions of the growth of fluctuations have been proposed: in it was 
suggested that: 

V L ~ C ■ V (5) 

with v ~ 0.5; on the other hand in (see also ||), a formula based on the dynamical 
renormalization calculation was proposed: 

Vz^C-LVMog^L/A)), (6) 
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with (3 = 1/4 and Lq = 8. In this section we will argue that this data supports the conclusion 
of [0], although we disagree on the identification of the constant L . 

Clearly it is difficult to distinguish between the behavior ([5]) and (^). For example, the 
quantity Vl/(L 1 / 2 logL) is a monotonic decreasing function of L throughout the range of our 
simulations and Vl/ (L 1 / 2 log 1 / 8 L) is monotonic increasing for L > 64, suggesting asymptotic 
behavior @ for some intermediate value of /3, but in fact the same is true if Vl is replaced 
by L 0-53 , the behavior found in ||. 

As a qualitative criterion for comparison of fit, we ask for what minimal value 2 K of the 
system size various proposed forms (all having two free parameters) can provide a good fit 
over interval [2 , 2 19 ], rejecting a fit as bad if it fails the standard x 2 test at a 99% confidence 
level. The asymptotic form fl5|), computed as a linear fit logtV^/L 1 / 2 ) ~ a + MogL, yields 
K = 12, while the forms (|6|), again as linear fits (Vl/ L 1 ^ 2 ) 1 ^ ~ ag + frglogL, yield K = 11, 
6, 8, and 10 for /3 = 1, 1/2, 1/4, 1/6, respectively. This analysis thus provides evidence to 
prefer (||), which describes the data over a wider range of system sizes, although it does not 
conclusively rule out (|). 

We now consider (|5|) and ask whether our data can determine the exponent (3. The 
difficulty is that the three parameters in (^) provide a great deal of freedom, allowing for 
good fits for a range of values of /3, if we vary the values of L and C. The test discussed 
above is indefinite, and in fact mildly favors a value of 1/2, as opposed to the renormalization 
group calculation of 1/4. We now turn to a different line of argument which does provide 
evidence for the latter value, and begin with a review of the second approximate approach 
of [|. 

As mentioned in the introduction, the application of the renormalization group depends 
on consideration of our dynamical model on the doubly infinite line, viewed as a height 
interface growth model. Such interface growth models may be described by non-linear 
stochastic diffusion equations of the EW and KPZ type; for our model the symmetries in 
the problem lead to the equation 
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dh d h dh ( dh\ . . . 

— = — ^ + — h u 3 — + ^(x, t) + higher order terms, (7) 
dt dx l dx \dx J 

where rj is a stochastic noise term. The term in (0) with coefficient v 1 is usually eliminated 
by a Galilean transformation, but this is not possible in the semi-infinite system and in fact 
this term has important consequences, discussed below. If v 3 = then (|7|) is the Edward- 
Wilkinson growth model, which leads to interface growth (h(t) 2 ) ~ t^. The full equation (|7|) 
can be analyzed using dynamical renormalization techniques f7|; in this analysis the higher 
order terms are irrelevant and the cubic term is marginally irrelevant, giving 

(h{t) 2 )~Athog*{t/t ). (8) 

The parameter v i is positive, so that fluctuations in h(t) travel to the right. We can compute 
the velocity v of infinitesimal density perturbations in the (Bernoulli) steady state, via 
v = dJ I 'dp\ p=1 , 2 , where J(p) is the current at density p, to find v = 8 j7|. To apply these 
observations to the original model on the half open system it was suggested in || that 
excitations created at the closed end of the system grow according to (^) as they travel 
towards the open end at speed v ; this leads to (^) with (5 = 1/4, C = A/v 1 ^ 2 , and L = vto. 
(See also @, although there it is tacitly assumed that to = !)• 

In order to understand the exact form of the logarithmic corrections and check the 
validity of the above prediction, we studied a class of modifications of the original model, 
with dynamics defined as follows: we choose a site at rate 1, then exchange the spin at that 
site with the k th spin of opposite sign to the right. The original model is obtained by taking 
k = 1. These models all share the symmetries of the original, so we expect that they will 
belong to the same universality class and hence have similar behavior of fluctuations; thus 
the variance V^k) in these models should satisfy (|6|) with the same value of (3 but with 
possibly different constants L (k) and C(k). We obtained data for models with k = 2, 3 
(see Table I). Again, we can fit (^) to the individual data sets for a range of values of (3 
if we choose appropriate L (k) and C(k); to proceed further, we would like some a priori 
argument determining how these values depend on k. 
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Now we propose a heuristic argument which describes how the constant C(k) might vary 
with k. On the one hand, a computation in the steady state determines that the current 
in this model is proportional to k and hence that the velocity v(k) of excitations satisfies 
v(k) = k ■ v(l) = 8k. On the other hand, the model for parameter k can be thought 
of as a modification of the original model in which, during one time step, k exchanges 
rather than one spin exchange take place. Now these k exchanges take place in correlated 
positions. These correlations are certainly significant on the microscopic scale (e.g., the 
distance between two successive moves would definitely be affected by these correlations), 
but as seen from the computation of the current and the velocity these correlations are 
apparently not important on the hydrodynamic scale. Hence there should be some length 
scale below which the correlations are significant. We have two parameters in (|6|) which 
could be affected by these correlations, C{k) and L (k). C{k) (along with (3) determines 
the behavior to leading order at large L, and Lo(k) is significant only at higher orders. If 
we assume that the correlations are not important at least to the order determining C(k), 
then for a computation of C(k), the model for parameter k would correspond to the original 
model with the time rescaled by a factor of k. This would imply that A(k) ~ A; 1 / 2 A(1), and 
hence that C(k) = A(k)/v(k) l l 2 = C(l). 

Another way of saying this, which does not make any reference to the time dependent 
problem, is to write (M£ +1 ) = (M|) + 2{a L+1 M L ) + 1. Since {M 2 L+l ) - (Ml) -» as 
L — > oo, we must have ((7l+iMl) — > — | + o(L). Now if the leading order term in o(L) 
is also independent of k (which seems not unreasonable) then we would also have C(k) 
independent of k and vice versa. 

To study the /^-dependence of C(k), we plot in Figure 1 (Vl(/c) / X 1 / 2 ) 1 /^ versus logL for 
(3 = 1/2, 1/4, and 1/6; according to equation (|6]) this curve is approaching a straight line 
with slope C(k) 1 ^ and offset C(k) 1 ^ log(L ). Only for [3 = 1/4 can a data collapse be 
obtained by vertically shifting the curves (which corresponds to a change in Lq), indicating 
that ^-independence of C{k) holds only for (3 = 1/4. Results of least-squares fit for all data 
on the interval [2 10 ,2 18 ], in the form (Vl/L 1 ' 2 ) 1 '^ ~ dp + feglogL, are given in Table 2; the 
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approximate constancy of &i/4(A;) verifies the conclusion. Thus under the assumption made 
above, this provides independent justification of the conclusion f3 = 1/4. 

We also observe from Table 2 that ap(k) and hence Lo(k) depends strongly on k for all 
(3. We interpret this as indicating that the correlations are important on the scale of Lq and 
hence we cannot use the above argument to determine the k dependence of L (k). 

To summarize, we now have numerical evidence for logarithmic corrections, and our data, 
when supplemented by a heuristic argument, supports the prediction (f|) (with (3 = 1/4) 
obtained with the help of dynamical renormalization techniques f7|. 



III. THE BIASED CASE 

In the biased case we generated data for L < 2 18 for A = 1/4. As compared to the data 
reported in ||, this provides much better evidence for the asymptotic form Vj, ~ BL 2 ^ 3 
predicted by both the KPZ and CVA approximations. We note that for the biased model 
the differential equation for the height evolution is the usual KPZ equation 

dh d 2 h dh ( dh\ 2 

— = v—— + Vi- — h t>2 — + T]{x, t) + higher order terms, (9) 
dt dx 2 dx \dx ) 

where now the the non-linear quadratic term is relevant and it changes the power of the 
growth of fluctuations in time from 1/2 to 2/3. 

We now compare the simulation results more closely with those obtained analytically for 
the asymptotic case of the CVA (the continuum limit): B cv = 1.544A 1/3 (1 - v / A) 2/3 (l + 
V\y 2 0. Note that B cv vanishes at A = 0, 1 and has a maximum at A = 7-4^ w 0.0718. 
As shown in Figure 2, the values obtained from the CVA seem to be strikingly close to 
the simulation values in the asymptotic limit. The data suggests that the constants of 
proportionality B for the variance in the two cases are extremely close if not identical. To 
check whether the CVA was close only for this particular value of A, we studied two more 
values of bias A = 1/2, 1/8. As shown in Figure 2 for both the cases, we observed that the 
CVA seemed to give values which are extremely close to the actual values. 
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In order to have a better estimate of the accuracy of the CVA, we also studied the whole 
distribution of magnetization. The continuum limit CVA gives for the distribution of the 
magnetization values the fourth power of an Airy function ||, with a somewhat ad hoc 
cutoff suggested by the discrete CVA (the simulations and the CVA both give a Gaussian 
distribution for A = 1). This function has an asymmetry about the mean value. In Figure 3 
we compare this CVA limiting distribution with the distribution of magnetization determined 
from simulations, with all distributions normalized to have a mean of and a standard 
deviation of 1. It can be seen that the actual distribution shows the same characteristic 
features as the CVA result; moreover, looking at the distributions for different L values 
suggests a slow approach of the distribution towards the Airy function with increasing size L. 
However, if we look at the leading measure of asymmetry in the distribution — the normalized 
third moment ^^^^2^3/2 — then we observe (see Figure 3 inset) that the actual values 
and that obtained from the CVA are quite different. 

We also studied the values of average magnetization obtained for finite size systems. As 
noted in |J the exact asymptotic value and that obtained from the CVA coincide and are 
given by (Ml) = /iL, where /z = In Figure 4, we plot the approach of the average 

magnetization to the asymptotic value. Both the approach for the CVA and that for the 
actual data show very similar behavior. 

At the end of this analysis, we are left with the somewhat puzzling result, that though 
the CVA does not reproduce, even qualitatively, the results for the unbiased case, it seems 
to be an extremely accurate description for the biased case. 
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FIGURES 

FIG. 1. If F L « C{k)L 1 / 2 log/ 3 (L/Lo), then (F L /L 1 / 2 ) 1 /' 3 as a function of log L is a straight line 
with slope C(k) 1 ^. In the figure, (V L /L 1 / 2 ) 1 / /3 is plotted as a function of log L for k = 1 (plusses) , 
A; = 2 (circles), and k = 3 (diamonds), and /3 = 1/4. The curves for At = 1 and A: = 2 are shifted 
vertically, to obtain a collapse. The insets show that for (3 = 1/2 (left inset) and (3=1/6 (right 
inset) the data does not collapse. 

FIG. 2. Vl/L 2 / 3 is plotted as a function of log L, for A = 1/8 (diamonds), A = 1/4 (circles), and 
A = 1/2 (squares). In the first two cases the data shows convergence to a constant, indicating that 
Vl ~ L 2 / 3 without logarithmic corrections. The lines are the CVA values for the corresponding 
biases; the asymptotic CVA values are 0.3150, 0.2723 and 0.1855, respectively. 

FIG. 3. Normalized distribution (with mean and standard deviation 1) of the magnetization 
Ml for the biased case (A = 1/4) with system sizes L = 1024 (dashed line) and L = 262144 
(dotted line) together with the normalized distribution obtained from the asymptotic CVA (solid 
line). The inset shows the third moment of the magnetization, as a function of logL, for the CVA 
(solid line) and for the simulation(circles). 

FIG. 4. Finite-size effects in the magnetization for A = 1/4: the deviation from (Ml) /L = 1/3 
versus system size L is given in a log-log plot. The solid line is the CVA, the circles are simulation 
results. 
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TABLES 



L 


V L = V L {1) 


Vl(2) 


Vl(3) 


2 


1.335[1] 


2 


2 


4 


1.918[1] 


2.75[1] 


3.637[1] 


8 


2.733[1] 


3.818[1] 


4.626[2] 


16 


3.890[1] 


5.271[1] 


6.451 [2] 


32 


5.576[2] 


7.338[3] 


8.868[6] 


64 


8.02[2] 


10.286[4] 


12.292[6] 


128 


11.58[3] 


14. 543 [4] 


17.20[1] 


256 


16.76[3] 


20.66[1] 


24.27[2] 


512 


24.22[4] 


29.44[2] 


34.37[2] 


1024 


35.01[3] 


41.96[6] 


48. 70 [3] 


2048 


50. 54 [9] 


59.94[8] 


69. 20 [3] 


4096 


73.1[1] 


85.6[2] 


98.41 [4] 


8192 


105. 3 [2] 


122.4[1] 


139.90[7] 


16384 


151. 7 [4] 


174. 6 [3] 


199 2 [21 


32768 


218.3[5] 


249. 9 [4] 


283.1 [3] 


65536 


313[1] 


355.4[4] 


402.4[5] 


131072 


449 [2] 


507[3] 


573 [2] 


262144 


644 [4] 


720 [2] 


813[3] 


524288 


924 [3] 




1162[5] 



TABLE I. The variances obtained for the unbiased case for different values of k. The quantity 
in brackets after each value is the statistical error in the least significant digit. 
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k 


a l/2 


°l/4 


a l/6 


1 


0.694 ± 0.018 


0.086 ± 0.050 


-0.969 ±0.102 


2 


1.376 ±0.023 


1.682 ± 0.084 


1.551 ±0.228 


3 


2.052 ± 0.018 


4.104 ± 0.084 


7 924 ± 300 


k 




h/4 


h/6 


1 


0.073 ± 0.002 


0.194 ±0.006 


0.386 ±0.013 


2 


0.050 ± 0.003 


0.184 ±0.009 


0.507 ±0.025 


3 


0.038 ± 0.002 


0.179 ±0.010 


0.640 ± 0.036 



TABLE II. Values of ap and bp obtained from least-squares fits on the data for the unbiased 
model. 
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L 


v L 


(M)l 


((M- < M >f) L 


v L 






* = * 


* = * 


A = i 


\ = \ 


^ = 1 


1024 


36.82[6] 


344.57[2] 


-0.171[6] 


35.7[1] 


37.0[2] 


2048 


55.3[1] 


686.98[2] 


-0.185[2] 


52.4[2] 


57.2 [3] 


4096 


83.4[2] 


1371.06[3] 


-0.204[4] 


76.4 [4] 


88.5 [6] 


8192 


128.8[4] 

L J 


2738.32[5l 


-0.225[3] 


111.9[5] 


136.6[8l 

L J 


16384 


193[2] 


5471.4[1] 


-0.240[3] 


165.5[6] 


215[2] 


32768 


302.7[6] 


10936.1[3] 


-0.256[2] 


247[1] 


344 [3] 


65536 


473 [1] 


21862.62[6] 


-0.266[4] 


370 [5] 


533 [5] 


131072 


735 [4] 


43713.2[2] 


-0.281[6] 


566 [5] 


849 [10] 


262144 


1165 [7] 


84105.9[2] 


-0.266[5] 







TABLE III. The variances for three different values of bias, and the first and third moments for 
A = j. Note that we have not obtained data for systems of size smaller than 1024. The statistical 
error in the least significant digit is given by the quantity in brackets. 
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